% Two type capital model

clear all
close all

global par g_T2
%========================================================================
% Initial variable value at T = 0:

initial.y0 = 1;     % output of goods - choose units so that it is normalized to 1.
initial.k0 = 3.2349;   % 3.6071282734*0.8967995, capital stock in goods and services production sector 
initial.kR0 = 0.3723;  % capital stock in the fossil fuel energy sector

initial.c0 = 0.6619974;       % consumption 
initial.g0 = 0.11068716;     % per unit mining cost
initial.n0 = 0.00832913;      % mining investment

%=========================================================================
% The constant parameters of the model:

%%% Lifetime utility function 

par.beta = 0.05;  % time discount factor
par.gamma = 4;    % the coefficient of relative risk aversion

%%% Production function y = Ak

par.A = 1/initial.k0;  

%%% Energy efficiency function

par.F0 = par.A;   % initial energy efficiency
par.Fbar = 0.05; % technological energy efficiency limit
par.a = 0.001;  % the slope 

%%% Differential equations of capitals

par.delta = 0.04; % the capital depreciation rate

%%% Energy production function 

par.muR = 1;  % fossil fuel energy production coefficient

%%% Mining cost function

par.Sbar = 2126.0527;           % Total feasible resources
initial.Res0 = 1./0.065100642;      % Initial estimate of proved reserves
initial.gS0 = 0.00015;              % Partial derivative of g with respect to S 

par.alpha3 = 15;               
par.alpha2 = (par.Sbar - initial.Res0).*par.alpha3;
par.alpha1 = initial.gS0.*initial.Res0.^2;
par.alpha0 = initial.g0 - par.alpha1./initial.Res0; 

%%% renewable energy production coefficient function
par.alpha = 0.25;
par.Gamma1 = (4*initial.g0).^(-1./par.alpha);
par.Gamma2 =0.8*initial.g0; 
par.m = 0.1163;  % operations and maintenance costs, calculated by data
par.H0 = 0 ;
par.Hbar = par.A*par.Gamma2/par.m;
par.b = 0.01;

%%% knowledge accumulation function
par.psi = 0.33;

%%% population function

par.Q0 = 1;
par.popgr = 0.01;

par.Abar = -par.beta-par.delta-(par.m*par.Fbar-par.A*par.Hbar)/(par.Hbar+par.Fbar);
%=========================================================================
%%% Make the guesses of unknown variables
T3 = 500;
k_T3 = 8000;

S_T2 = 2000;
N_T2 = par.alpha2/(par.Sbar-S_T2)-par.alpha3+0.5;
kR_T2 = 500;

%%% Regime 4: Analytical solution

kB_T3 = par.Fbar*k_T3/par.Hbar;
lambda_T3 = ((par.Hbar+par.Fbar)*(par.beta*par.gamma+par.Abar* ...
                                  (par.gamma-1))*k_T3/(par.gamma*par.Hbar))^(-par.gamma);


%%% Regime 3:  Renewable only, iR=0, n=0, j > 0, iF >= 0

h_T3 = (par.Hbar - par.H0) / par.b;
f_T3 = (par.F0 - par.Fbar) / par.a;

eta_T3 = 0;
phi_T3 = 0;

% $$$ events3 g_T2 = (A-delta)/(ak+F0-af) detects that rho changes
% $$$ from 0 to 1. Fossil fuel is used after this point.
g_T2 = par.alpha0 + par.alpha1/(par.Sbar-par.alpha2/(par.alpha3+N_T2)-S_T2);

z_T3 = [k_T3, lambda_T3, kB_T3, h_T3, f_T3, phi_T3, eta_T3];
tspan3 = [T3 0];

% options = odeset('RelTol',5e-14,'AbsTol',5e-14,'events',@events2);
% [t3,z3,TE2,ZE2,IE2] = ode45(@diff_Reg3, tspan3, z_T3, options);
% $$$ For test purpose, solve the ODE without the events:
[t3,z3] = ode45(@diff_Reg3, tspan3, z_T3);


% $$$ Extract components of z for graphing

k3 = z3(:,1);
lambda3 = z3(:,2);
kB3 = z3(:,3);
h3 = z3(:,4);
f3 = z3(:,5);
phi3 = z3(:,6);
eta3 = z3(:,7);

% % $$$ Extract initial values to use as terminal values in the following
% % $$$ regime. 
% 
% T2 = TE2;
% k_T2 = ZE2(1);
% lambda_T2 = ZE2(2);
% kB_T2 = ZE2(3);
% h_T2 = ZE2(4);
% f_T2 = ZE2(5);
% % $$$ phi_T2 = ZE2(6);
% eta_T2 = ZE2(7);
% 
% 
% % % %%% Regime 2: Fossil fuel use phases out, iR>=0, n>=0, j>0, iF>0.
% % 
% % sigma_T2 = 0;
% % nu_T2 = 0;
% % qR_T2 = 0;
% % 
% % z_T2 = [k_T2, lambda_T2, kB_T2, h_T2, f_T2, eta_T2, kR_T2, S_T2, N_T2, ...
% %         sigma_T2 nu_T2, qR_T2];
% % tspan2 = [0 T2];
% % 
% % options = odeset('RelTol',5e-14,'AbsTol',5e-14,'events',@events1);
% % [t2,z2,TE1,ZE1,IE1] = ode45(@diff_Reg2, tspan2, z_T2, options);
% % 
% % %%% Test if h = 0.
% % 
% % 
% % % $$$ Extract components of z for graphing
% % 
% % k2 = z2(:,1);
% % lambda2 = z2(:,2);
% % kB2 = z2(:,3);
% % h2 = z2(:,4);
% % f2 = z2(:,5);
% % eta2 = z2(:,6);
% % kR2 = z2(:,7);
% % S2 = z2(:,8);
% % N2 = z2(:,9);
% % sigma2 = z2(:,10);
% % nu2 = z2(:,11);
% % qR2 = z2(:,12);
% % 
% % % $$$ Extract initial values to use as terminal values in the following
% % % $$$ regime. 
% % 
% % T1 = TE1;
% % k_T1 = ZE1(1);
% % lambda_T1 = ZE1(2);
% % f_T1 = ZE1(5);
% % kR_T1 = z(7);
% % S_T1 = z(8);
% % N_T1 = z(9);
% % sigma_T1 = z(10);
% % 
% % %%% Regime 1: Fossil fuel only. n>0, iR>0.
% % 
% % z_T1 = [k_T1, lambda_T1, f_T1, kR_T1, S_T1, N_T1, sigma_T1];
% % tspan1 = [0 T1];
% % 
% % options = odeset('RelTol',5e-14,'AbsTol',5e-14);
% % [t1,z1] = ode45(@diff_Reg1, tspan1, z_T1, options);
% % 
% % 
% % 
% % 
% % % $$$ Extract components of z for graphing
% % 
% % k1 = z1(:,1);
% % lambda1 = z1(:,2);
% % f1 = z1(:,3);
% % kR1 = z1(:,4);
% % S1 = z1(:,5);
% % N1 = z1(:,6);
% % sigma1 = z1(:,7);
% % 
% % k0_diff = k1(length(k1)) - initial.k0;
% % kR0_diff = kR1(length(kR1)) - initial.kR0;
% % S0 = S1(length(S1));
% % N0 = N1(length(N1));
